# Manuscript Replication File
# Trojan Horses in Liberal International Organizations?:
# How Democratic Backsliders Undermine the UNHRC
# September 2023

# Load packages ####
library(dplyr)
library(ggplot2)
library(scales)
library(tibble)
library(nnet)
library(AER)
library(broom)
library(MASS)
library(sandwich)
library(lmtest)
library(topicmodels)
library(readtext)
library(stm)

# Figure 1: Average electoral democracy index scores of UNHRC members by year ####
rm(list=ls())
load("unhrc_members_data.RData")

unhrc_members %>%
  group_by(year) %>%
  summarise(mean_polyarchy = mean(v2x_polyarchy, na.rm = T)) %>%
  ungroup() %>%
  filter(year > 1989) %>%
  ggplot(aes(y = mean_polyarchy, x = year)) + 
  geom_smooth(col = "steelblue4", se = FALSE) + 
  theme_bw() + 
  xlab("") + 
  ylab("Average Electoral Democracy Index") + 
  theme(text = element_text(size=20))

# Figure 2: Distribution of Votes in UNHRC, 2006–2021 ####
rm(list=ls())
load("unhrc_regression_data.RData")

# Figure 2a: 2006-2011
unhrc_regress %>%
  filter(vote != 4 & year >=2006 & year < 2012) %>%
  ggplot(aes(x = vote)) +
  geom_bar(aes(y = (..count..)/sum(..count..)), fill = "steelblue4") + 
  theme_bw() + 
  xlab("") + 
  scale_x_discrete(labels=c("1" = "No", "2" = "Abstain",
                            "3" = "Yes")) + 
  labs(title = "", 
       x = "", y = "") + 
  scale_y_continuous(labels=percent, limits = c(0, 1))  +
  theme(text = element_text(size=20))

# Figure 2b: 2012-2021
unhrc_regress %>%
  filter(vote != 4 & year >= 2012) %>%
  ggplot(aes(x = vote)) +
  geom_bar(aes(y = (..count..)/sum(..count..)), fill = "steelblue4") + 
  theme_bw() + 
  xlab("") + 
  scale_x_discrete(labels=c("1" = "No", "2" = "Abstain",
                            "3" = "Yes")) + 
  labs(title = "", 
       x = "", y = "") + 
  scale_y_continuous(labels=percent, limits = c(0, 1)) +
  theme(text = element_text(size=20))

# Figure 3: Percent of yes votes in the UNHRC by year, 2006–2021 ####
rm(list=ls())
load("unhrc_regression_data.RData")

yes_votes <- unhrc_regress %>%
  group_by(year) %>% 
  count(vote, keep.all = T) %>%
  mutate(total = sum(n)) %>%
  mutate(prop = (n/total)) %>%
  filter(vote == 3) %>%
  ungroup()

ggplot(data = yes_votes) + 
  geom_line(aes(x = year, y = prop), col = "steelblue4") +
  labs(title = "", 
       x = "", y = "") + 
  theme_bw() + 
  theme(text = element_text(size=20)) +
  scale_x_continuous(breaks=c(2006,2011,2016, 2021),
                   labels=c(2006,2011,2016, 2021)) +
  scale_y_continuous(labels=percent)

# Figure 4: Total number of votes in the UNHRC per country, 2006–2021 ####
rm(list=ls())
load("unhrc_regression_data.RData")

# calculate number of votes by type per country
stats <- subset(unhrc_regress, vote != 4)
tab <- table(stats$votingstate, stats$vote)
tab <- as.data.frame.matrix(tab)
tab <- tibble::rownames_to_column(tab, "Country")
colnames(tab)[2] <- "no"
colnames(tab)[3] <- "abstain"
colnames(tab)[4] <- "yes"

# table of total votes per country
tab <- tab %>%
  mutate(total = no + abstain + yes) %>%
  dplyr::select(Country, total)

# plot
tab %>%
  ggplot(aes(y = Country, x = total)) + 
  geom_bar(stat = "identity", color = "gray", fill = "steelblue4") + 
  theme_bw() + 
  xlab("") + 
  ylab("") + 
  ggtitle("")

# Table 1: Democratic backsliding and vote choice in the UNHRC, 2006–2021 ####
rm(list=ls())
load("unhrc_regression_data.RData")

# drop obs where vote == 4 (missing)
mod <- subset(unhrc_regress, vote != 4)

# set baseline level of dv ('yes' vote) and relevel
mod$vote <- ifelse(mod$vote == "1", "no", mod$vote)
mod$vote <- ifelse(mod$vote == "2", "abstain", mod$vote)
mod$vote <- ifelse(mod$vote == "3", "yes", mod$vote)
mod$dv <- as.factor(mod$vote)

# reference category 
mod$dv <- relevel(mod$dv, ref = "yes")

# multivariate logit with year fixed effects
tab1.mod1 <- multinom(dv ~ db_bin_democracies + 
                                               v2x_polyarchy_votingstate +
                                               Israel + 
                                               Syria +
                                               v2x_polyarchy_targetstate + 
                                               same_geo_region + 
                                               vid_usagree + 
                                               tid_usagree + 
                                               vid_tid_unga_agree + 
                                               tid_hr_score +
                                               factor(year), data = mod)

# clustered standard errors 
se.tab1.mod1 <- coeftest(tab1.mod1, vcov = vcov(tab1.mod1, type = "HC0", cluster = mod$vid))

# relative risk ratios
exp(coef(tab1.mod1))

# Figure 5: The percent of UPR reviewers that qualify as cases of backsliding has steadily increased over time ####
rm(list=ls())
load("upr_regression_data.RData")

upr_regress %>%
  group_by(ccode_recommending_state, year) %>%
  filter(row_number() == 1) %>%
  group_by(year) %>%
  summarise(num_db = sum(db_bin_democracies), 
            num_reviewers = n_distinct(ccode_recommending_state), 
            percent_db = (num_db/num_reviewers)*100) %>%
  ungroup() %>%
  ggplot(aes(x = year, y = percent_db)) + 
  geom_bar(stat = "identity", fill = "steelblue4") + 
  theme_bw() + 
  xlab("") + 
  ylab("Percent of Backsliders among UPR Reviewers") + 
  theme(text = element_text(size=20))

# Table 2: UPR Issues Identified against Advanced Democracies, 2008–2020 ####
rm(list=ls())
load("upr_regression_data.RData")

# subset to only include obs where advanced democracy is under review & reviewer is not adv demo
mod <- subset(upr_regress, upr_regress$advanced_demo_state_under_review == 1 & 
                   upr_regress$advanced_demo_rec_state == 0)

# model 1: all reports 
tab2.mod1 <- glm(dv_count_issues_per_review ~ db_bin_democracies +
                   v2x_polyarchy_recommending_state + 
                   v2x_polyarchy_state_under_review + 
                   same_geo_region + 
                   usagree_recommending_state + 
                   usagree_state_under_review + 
                   review_recommend_unga_agree +
                   hr_score_state_under_review + 
                   recommending_under_review + 
                   factor(year), family = "quasipoisson", data = mod)

# model 1 clustered standard errors
(se.tab2.mod1 <- coeftest(tab2.mod1, vcov = vcovHC(tab2.mod1,type = "HC0", 
                                                   cluster = mod$ccode_recommending_state)))

# model 1 dispersion parameter 
dfun <- function(object) {
  with(object,sum((weights * residuals^2)[weights > 0])/df.residual)
}
dfun(tab2.mod1)

# model 2: reports by UNHRC members

# subset to only include observations where reviewer is UNHRC member
mod2 <- subset(mod, mod$unhrc_member_recommending_state == 1)

tab2.mod2 <- glm(dv_count_issues_per_review ~ db_bin_democracies +
                   v2x_polyarchy_recommending_state + 
                   v2x_polyarchy_state_under_review + 
                   same_geo_region + 
                   usagree_recommending_state + 
                   usagree_state_under_review + 
                   review_recommend_unga_agree +
                   hr_score_state_under_review + 
                   recommending_under_review +
                   factor(year), family = "quasipoisson", data = mod2)

# model 2 clustered standard errors
(se.tab2.mod2 <- coeftest(tab2.mod2, vcov = vcovHC(tab2.mod2, type = "HC0",
                                                   cluster = mod2$ccode_recommending_state)))

# model 2 dispersion parameter 
dfun(tab2.mod2)

# UPR STM Data Prep and Models ####
rm(list=ls())
load("upr_stm_data.RData")

# Cleaning and Preprocessing ####
# drop observations where any covariates are missing
upr_stm <- subset(upr_stm, !is.na(upr_stm$db_bin_democracies))
upr_stm <- subset(upr_stm, !is.na(upr_stm$ccode_recommending_state))
upr_stm <- subset(upr_stm, !is.na(upr_stm$v2x_polyarchy_recommending_state))
upr_stm <- subset(upr_stm, !is.na(upr_stm$v2x_polyarchy_state_under_review))
upr_stm <- subset(upr_stm, !is.na(upr_stm$same_geo_region))
upr_stm <- subset(upr_stm, !is.na(upr_stm$usagree_recommending_state))
upr_stm <- subset(upr_stm, !is.na(upr_stm$usagree_state_under_review))
upr_stm <- subset(upr_stm, !is.na(upr_stm$hr_score_state_under_review))
upr_stm <- subset(upr_stm, !is.na(upr_stm$recommending_under_review))
upr_stm <- subset(upr_stm, !is.na(upr_stm$review_recommend_unga_agree))

# create corpus
processed <- textProcessor(upr_stm$recommendations, metadata = upr_stm)
out <- prepDocuments(processed$documents, processed$vocab, processed$meta)
docs <- out$documents
vocab <- out$vocab
meta <- out$meta

# Estimate STM Model with 6 topics ####
stm <- stm(docs, vocab, K = 6, 
                    prevalence = ~ db_bin_democracies + 
                      v2x_polyarchy_recommending_state + 
                      v2x_polyarchy_state_under_review + 
                      same_geo_region + 
                      usagree_recommending_state + 
                      usagree_state_under_review + 
                      review_recommend_unga_agree +
                      hr_score_state_under_review + 
                      recommending_under_review + factor(year),
                      max.em.its=200, data=out$meta, init.type="Spectral", 
                      seed=84581)

# Figure 6: Labels for STM Topics ####
labelTopics(stm)

# Figure 7: STM Model Coefficient Plot (binary measure of backsliding) ####
# estimate effects
stm.effect <-  estimateEffect(1:6 ~ db_bin_democracies +  
                                          v2x_polyarchy_recommending_state + 
                                          v2x_polyarchy_state_under_review + 
                                          same_geo_region + 
                                          usagree_recommending_state + 
                                          usagree_state_under_review + 
                                          review_recommend_unga_agree +
                                          hr_score_state_under_review + 
                                          recommending_under_review +
                                          factor(year), stm, meta=processed$meta, uncertainty = "Global")

# plot effects for binary measure of backsliding 
par(mar=c(7,9,7,1)+.1)
plot(stm.effect, 
     covariate = "db_bin_democracies", 
     topics = c(1:6), 
     method = "difference", 
     cov.value1=1, # set the treatment to cov.value1
     cov.value2=0, # set as the comparison group
     ci.level = 0.95, 
     xlab = "Change in topical prevalence from non-backsliding (0) to backsliding (1)",
     main = "", labeltype = "custom", 
     custom.labels = c("Human development",
                       "Ratify international human rights conventions",
                       "Women and LGBTQ+ rights", 
                       "Civil and political liberties",
                       "Establish national human rights institutions",
                       "Protect minorities and vulnerable groups"),
     cex.lab=1.5, cex.axis=1.5)

# Figure 8: STM Model Coefficient Plot (continuous measure of regime type) ####
par(mar=c(7,9,7,1)+.1)
plot(stm.effect, 
     covariate = "v2x_polyarchy_recommending_state", 
     topics = c(1:6), 
     method = "difference", 
     cov.value1=0.33, 
     cov.value2=0.87,  
     ci.level = 0.95, 
     xlab = "Change in topical prevalence when EDI declines (3rd to 1st quartile)",
     main = "", labeltype = "custom", 
     custom.labels = c("Human development",
                       "Ratify international human rights conventions",
                       "Women and LGBTQ+ rights", 
                       "Civil and political liberties",
                       "Establish national human rights institutions",
                       "Protect minorities and vulnerable groups"),
     cex.lab=1.5, cex.axis=1.5)


